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ABSTRACT 

A numerical model of the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect for ob- 
jects defined in terms of a triangular mesh is described. The algorithm requires that each 
surface triangle can be handled independently, which implies the use of a ID thermal model. 
Insolation of each triangle is determined by an optimized ray-triangle intersection search. Sur- 
face temperature is modeled with a spectral approach; imposing a quasi-periodic solution we 
replace heat conduction equation by the Helmholtz equation. Nonlinear boundary conditions 
are handled by an iterative, FFT based solver. The results resolve the question of the YORP 
effect in rotation rate independence on conductivity within the nonlinear ID thermal model 
regardless of the accuracy issues and homogeneity assumptions. A seasonal YORP effect in 
attitude is revealed for objects moving on elliptic orbits when a nonlinear thermal model is 
used. 

Key words: radiation mechanisms: thermal — methods: numerical — celestial mechanics — 
minor planets, asteroids 



1 INTRODUCTION 

Referring to p art ial results of his predecessors (most notably 



Keternne to p art ial results or nis predecessors (most notably 
|Paddackril969l) ). iRubincamI ilOOd) forged the acronym 'YORP 
effect' (Yarkovsky-O'Keefe-Radzievskii-Paddack) to describe the 
influence of radiation effects on the rotation of a Sun-orbiting ob- 
ject. The radiation incident on the surface of a celestial body acts 
in three different ways: by the direct pressure, by the recoil force of 

reflected photons, and by the thermal radiation force. 

According to a simple geometric argument of 'Rubincam' 

ilOOd). further elaborated by iNesvornv & Vokrouh lickv (2008b) 



where Aii is the scattered fraction of incident power flux ~ the 
product of albedo A and of the flux E hitting the infinitesimal sur- 
face element d5 , and 



M, = -^^etr4rxdS, 



(3) 



Rubincam & Paddackl Jioi oT the average torque due to direct 



where the re-radiated energy flux is the product of the Stefan- 
Boltzmann constant cr, grey body emissivity factor ft and the fourth 
power of surface temperature T. Both terms are divided by the ve- 
locity of light c. 

Conservation of energy on the body surface implies that 



and 

radiation pressure vanishes. Physical properties of asteroid surfaces 
do not suggest a significant contribution of specular reflection, so 
we may focus on the remaining two phenomena: scattered (i.e. dif- 
fusively reflected) radiation and thermal re-radiation, defining the 
YORP torque M as the sum 



eto-r'* + a:m • vr - (1 -A)£ = 0, 



(4) 



M = M(j + Mu 



(1) 



of the torque Md generated by the scattered radiation, and of the 
grey body thermal radiation torque Mt. 

Within our present Lambertian model, the primary definitions 
of the YORP torque components in reference frame attached to an 
object's centre of mass are given as integrals over the body surface: 



so the absorbed flux {l-A)E is distributed between the re- 
radiation, proportional to T'*, and conduction term given as the 
product of thermal conductivity K and the normal derivative of tem- 
perature - the gradient projected on the outward normal unit vector 
n. 

Substituting the boundary condition l|4) into Eq. l[3]l we can 
merge a part of Mt with M^, so that the YORP torque becomes the 
sum 



M = Mr + Mc, 

of the principal term 



3 c 



AErxdS, 



(2) 
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Mr = -— ^^ErxdS, 

and the complement due to conductivity 



Mc= — (D K(n-VT)rx dS. 

3c Js 



(5) 



(6) 



(V) 



The subscript R refers to the usual 'Rubincam approximation' 
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of zero conductivity, and the problem of finding Mr is actually 
an exercise in computational geometry. Its most difficult part is 
the evaluation of E, discussed in Sect. [3] Computing the conduc- 
tivity term requires solving the heat diffusion equation. Our 
simplified ID thermal model is presented in Sect. |4] Two of its 
possible extensions are given in Appendices |B] and |C] but the lat- 
ter serves mostly as a theoretical argument and has not been im- 
plemented. The results of test runs with asteroids 1998 KY26 and 
6489 Golevka are presented in Sect.|5] In our opinion, they reveal a 
previously unnoticed seasonal YORP effect in attitude. Additional 
assumptions of our model are enumerated in Sect. |2l but we hope 
to relax them in future. 



2 PRELIMINARIES 
2.1 Body shape model 

Although there are many possible ways to describe the shape of a 
celestial body, the YORP studies practically rely on two variants: 
a spherical harmonics model (typical for analytical considerations) 
or a triangular mesh. We adopt the latter as more general, capable 
of representing even very irregular objects, and more suitable for 
the occlusion tests. As a consequence, an integral over the body 
surface becomes the sum of cubatures over all triangular patches 
forming the mesh. Of course, the real information about the surface 
points is given only at the vertices r,-, so the values of distance or 
any other coordinates dependent function have to be interpolated 
on a patch. In principle, it can be done using various interpolation 
rules, even the ones that involve the whole set of vertices, but in 
the YORP practice all authors rely on the local, linear interpolation, 
considering flat triangles and replacing all cubatures over triangular 
patches by the first order Gaussian midpoint rule 



e-, = b'= e. 



©Sun 



/(r)d5,«/(r,)Sj 



(8) 




perihelion 



equinox 



Figure 1. Reference frames and angles in a body-centered system. 



Thus, for a model of a celestial body with A'^i triangular faces, 
the YORP torque is approximated as a sum 



M 



with 



(13) 



(14) 



where two terms in the bracket are: Ej ~ the incident power flux 
evaluated at the centroid r j, and 



(15) 



These two terms are responsible for the Rubincam part and the con- 
ductivity complement, respectively. 



where S j is the area of a triangle determined by vertices r^, Tj, r^. 



and Yj is the centroid 



'■i = 7(l') + '-i+'-2)- 



(9) 



In particular, the oriented surface vectors S j = iij S j are constant on 
each triangular face, easily computed as 



(10) 



Of course, the mesh should be properly oriented, so that S j 
computed from Eq. l llOt is always directed along the outward nor- 
mal. The routine tests rely on checking the Gauss identity 



V5, = 0, 



(11) 



followed by asking if the volume resulting from the sum of oriented 
tetrahedral simplices 



(12) 



7=1 



is positive, when all A'm faces are included. Yet, even if both tests 
have been passed, there remains a number of possible degeneracies, 
like edges shared by more than two triangles, duplicated vertices 
etc., that are best to be checked before using the mesh. 



2.2 Dynamics 

Although recent works of IVokrouhlickv et il\ ( |2007|) and 
ICicalo & ScheeresI ( I2OI0I) have revealed the importance of tum- 
bling rotation for the dynamics under the YORP torque, we adhere 
to the usual assumption of the rotation around the principal axis of 
inertia - the e, unit vector of the body-frame basis. The principal 
axis mode remains a decent approximation o ver significant frag- 
ments of the evolutionary paths presented by IVokrouhlickv et aU 
and in some instances it may be possible to incorporate 
tumbling by rotating the basis with respect to the principal axes. 
So, we consider the principal axis mode equations of motion for 
the rotation rate oj, the obliquity e (the angle between the normal 
to the orbital plane and the spin vector o), parallel to e^), and the 
sidereal time Q (measured from the ascending node of the Sun on 
the object's equator to the body-frame basis e.^ vector) 



C ' 
M-ei 



Q. 



M-e2 



(16) 
(17) 
(18) 



ojC tane 

where C designates the maximum moment of inertia in the princi- 
pal axes frame. 

According to the above equations, the dynamics is governed 
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by the components of the YORP torque M in another kind of equa- 
torial reference frame (ej , 62, 63) (see Fig.lTJ with the same (centre 
of mass) origin as the body frame (e.v, e,., e^), with the same z di- 
rection, but with the re maining axes relate d to the equinox instead 
of to the principal axes ( lBreiteretall2010h 

e\ = smile x + co?,Q.ey, 

62 = -cosfle.v + sinflgj,, (19) 
^3 = 

2.3 Average YORP effect 

From the point of view of long term, systematic influence of the 
YORP effect, we are mostly interested in the mean values of a), 
s and Cl, with all daily and orbital periodic effects averaged out. 
Thus, assuming a uniform rotation with constant frequency oj and 
the Keplerian motion around the Sun, we need to find the mean 
values, defined for any function / 

</> = TT /d^dQ, (20) 



E = 0^(r, «, hQ)h- Mq, 



(25) 



47T^ Jo 

where £ is the mean anomaly of the Sun. We perform this averag- 
ing indirectly, working with discrete Fourier transforms (DFT) of 
f(C,Ci). Using the conventions described in Appendix [A] we as- 
sume 



</>=/[0], 



(21) 



where / = F2/ is the DFT of the function /, and / is the vector of 
samples of this function at different Cl and { values. 

According to Eq. il4i . the mean value of M ■ ^3 on a single 
triangle j is simply 



MJ-es) = -— {Ej[0]+Qj[0]) {rjxSj) e,, 



(22) 



where Sj is computed in the body frame, so the scalar product with 
63 = e- amounts to selecting the third component. But the averag- 
ing of the next two terms is slightly more involved: confronting 
Eq. ([19J with {AgJ we find 



- [%{Ej[1]+Qj[1]) {rjxSj)-ey 



-5{Ejll]+Qjll]) {rjxSj)-e, 



(23) 



(24) 



{Mj-e2) = -[%{Ej[l]+Qj[l]){rjxS,)-e., 
+5{Ej[l]+Qj[l]) {rjxSj)-e,]. 

The above expressions make use of the Hermitian property of the 
DFT for a real-valued function, and assume that in the sampling 
described in Appendix lAl the first angle is <^ = ^, and the second is 
i// = Sl. 



3 RUBINCAM TERMS 
3.1 Insolation function 

The major difficulty in dealing with the Rubincam part of the YORP 
effect is the computation of E, known as the insolation or irradia- 
tion function. In principle, the incident energy flux hitting a given 
surface element is a sum of two components: the direct flux from 
the Sun, and the radiation exchange complement, i.e. the energy 
coming from other elements of the body surface (either reflected 
or re-emitted). In the present paper we adhere to the approximation 
used in all previous works and consider only the direct part 



where O designates the solar radiation power flux at a given dis- 
tance of the body from the Sun r^. Using the solar constant Oo ~ 
1366 Wm"^, we have 



(26) 



with the reference distance (Iq = 1 au. 

Defining and computing the visibility function f is the heart 
of the problem. Its values are either f = 0, when the Sun is not vis- 
ible over the current surface element, or 1 otherwise. For convex 
bodies, ^ depends only on the scalar product of the outward normal 
unit vector h and the unit vector directed to the Sun Mq. In other 
words, whenever the zenith distance of the Sun is less than 90°, 
the visibility function f = 1, because the formal horizon (a local 
tangent plane perpendicular to h) and the actual horizon (the part 
of a celestial hemisphere not occluded by other surface elements) 
coincide for a convex object. In this case, computing ^ is so cheap 
and easy, that often the bodies of an arbitrary shape are formally 
treated as co nvex when com puting E, which is necessary in analyti- 
cal theories jNesvornv & Vo krouhlickv 200 ^2008j : lMvserJl2008l: 
iBreiter & MichalskalbOojTlBreit er et al. 201^ and handy in nu- 
merical or some semi-ana l ytical models ( Vokrouhlick v et al. 20o3; 
IScheeres & Gaskellll2008l: fcicalo & Scheereslboid) . However, the 
weakness of such pseudo-convex treatment for irregular, bouldered 
and cratered objects, test ified by Scheeres et al. (2008) and - in a 
quite different form - by Breiter et al. ( 2009), suggests to avoid it 
whenever possible, unless the shape model is known in advance to 
be convex (e.g. when it comes from the convex lightcurve inversion 
algorithm). 

Leaving aside the visibility function algorithms, to be dis- 
cussed in next subsections, we begin computations with tabulating 
the flux O and the components of Mq in the orbital frame for the 
mean anomaly i sampled at A' equidistant points in the full angle 
range < f < 27r. In the orbital frame, the direction cosines are for- 
mally 



OJorb ■ 



COs(Wo+/o) 

sin(a;o+/o) 




(27) 



where oJq is the argument of perihelion and /□ is the true anomaly 
of the Sun. Thus the two nonzero components and O are tabulated 
once, before the the main loop over surface triangles begins, so 
the cost of solving Kepler equation is relatively negligible. Other 
quantities precomputed before the main algorithm starts are: cen- 
troid positions rj, areas S j, and unit normal vectors nj associated 
with each triangular face. 

Given a pair of mean anomaly and rotation phase ({, Q.) we 
transform the solar vector to the body frame by means of two rota- 
tions: around the first axis by angle (-£), and then around the third 
axis by angle Cl, so that 



M0 



COS Cl COSE sin 
-sinfi cos£Cos£l 
sine 



where '*' are placeholders for irrelevant matrix entries. 



(28) 



3.2 Visibility function 

The fundamental operation in the evaluation of the visibility func- 
tion is the 'stabbing query', i.e. testing the intersection of a ray 
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ri>0. 



(29) 



with other surface triangles. T iiis standard tool of compu tational 
geometry is well documented jMoUer & Trumbord [19971) , so we 
skip the details focusing on a less trivial question: how to minimize 
the number of its calls. 

Obviously, there is no need to perform the query when the 
Sun is below the formal horizon, i.e. when nj • Wq < 0. So, the most 
straightforward selection tool is to create for each j a list of all 
triangles with at least one vertex above the formal horizon and, 
if ■ «Q > 0, perform the queries only with triangles from the list. 
But expecting that the list should be short is wishful thinking, based 
upon a false intuition of a flat landscape with distant mountains on 
the horizon and plenty of clear sky above a spectator's head. 

Quite a number of difficult to spot errors may arise if an opti- 
mized visibility algorithm is created with such a picture in mind. If 
there are craters or boulders on an asteroid, one should rather try to 
imagine the landscape seen by an ant climbing a pit or walking on 
a side of a boulder surrounded by other rough terrain features. The 
region of a clear sky can be a small, irregular, non-convex area, and 
its intersection with the daily Sun path can be a union of disjoint 
segmentsQ For a triangle on a boulder or a crater side, up to 90% 
of remaining triangles may stretch above the formal horizon, and it 
means that the number of queries should be additionally optimized. 

An o ptimization method, very briefly reported in the paper 
of IStatleri t2009.) . relies on a horizon map - a ID array of maxi- 
mum elevation of surface features above the formal horizon on a 
grid of local azimuth values for a given triangle. However, it is not 
clear from the author's description how far his approach is based 
upon the 'hills on a horizon' paradigm and whether he avoided the 
problems arising when some triangles overhang the local zenith, 
because then the altitude of clear sky has both the lower and the 
upper bound (smaller than 90°). 

Another, more robust way to hand le the op t imiza- 
tion, applied in papers like jv okrouhlic kv & CapeM I2OO2I : 
ICapek & Vokr ouhlick'(ll2004l : lDurech et aU200^ , and described in 
(Capek 2008, Appendix B2), amounts to creating a huge collection 
of 2D visibility tables for a specific object. For each surface 
triangle, a longitude-latitude Mercator map with 0/1 values on 
a 1° X 1° grid is first computed and stored in a file. During the 
YORP computation, the longitude and latitude of the Sun are 
rounded to full degrees and compared with the related entry of 
the visibility table. Creating the tables is time consuming, but 
performed only once for a given object. The drawbacks are: fixed 
discretisation error with uneven resolution on a sphere, and huge 
file space requirements. The largest shape model attacked with this 
approach was the t riangulated Itokawa shape with 196608 facets 
jfiurech et alll2008h . 

Our approach is an attempt to jo in the robustness of 
with the potential efficiency of IScheeresI ( |2007|) or 



Capek 



Statlei 



2009 '). The method has already proved its valor in computing 
the Rubincam part of the YORP effect for Itokawa and Eros us- 
ing t heir models of over 3 x 10^ triangular faces jBreiter et al.l 
20091) . Computing the YORP torques for a number of obliq- 
uity values £,-, we use the following arrangement of loops: tri- 
angles(obliquity(orbit(rotation)). The efficiency of our approach 



hinges upon the possibility of considering surface elements one by 
one in the outermost loop, which is possible within the assumptions 
of the illumination and thermal model of the present work, although 
suppressing the present restrictions in future, we will most likely 
find ourselves in a less comfortable situation. 

For a current surface triangle j we first create a 'horizon ar- 
ray', partitioning the local hemisphere into a fixed number of sec- 
tors (along meridians) and zones (along constant altitude circles). 
A typical setup uses about 100 sectors and 64 zones. Each triangle 
above the formal horizon is centrally projected onto the unit sphere 
with the origin at rj in order to find its 'bounding box' in azimuth 
and altitude. The problem has its own subtleties: the extreme az- 
imuth values are those of the vertices, but care must be taken about 
the cases of crossing the zero meridian; on the other hand, the ex- 
tremes of altitude are often different from the altitude of vertices 
due to the bending of a straight edge in central projection. And, last 
but not least, if some triangle is intersected by the local zenith line 
(parallel to the normal vector hj), it should be marked as a 'zenith 
triangle' and requires a special treatment, having a constant altitude 
circle as the bounding box. 

In addition to the bounding box determination, each triangle 
k is also labeled as foreground or background object, depending 
on the sign of the scalar product of its outward normal vector % 
and the relative position of its centre rj. -rj. Obviously, if any ray 
w intersects a foreground triangle, it must intersect a background 
triangle as well, hence ~ for economy of time and storage - we 
select a less populated of the foreground and background subsets of 
faces above the formal horizon as the candidates for future stabbing 
queries. 

Once the first loop over triangles kt j is completed, the hori- 
zon array is dynamically created with an appropriate size. Then, in 
the second loop over k + j, the number k is stored in the lists re- 
ferring to all zone-sector cells covered by the bounding box of the 
k-th triangle. 

Thus we create the horizon array - a set of lists containing 
possible occluders for a given solar altitude and azimuth. Actu- 
ally, the array covers the entire hemisphere only in the presence 
of a zenith triangle. If no such face is detected, we record the the 
bounding altitude of the clear sky cap and set the horizon array 
cells subdividing only the sector between the formal horizon and 
the clear sky limit circle. After that, the remaining computations 
are straightforward: fixing the value of obliquity (or opening the 
obliquity loop) we sample the rotation phase and mean anomaly, 
and for each pair of these angles compute the Sun vector Mq. If the 
Sun is above the formal horizon, we select an appropriate entry of 
the horizon array and perform stabbing queries with triangles from 
the list, until we record the intersection (f = 0) or the end of the 
list is reached (f = 1). Having collected all values of the insola- 
tion function E j e M.^~ , we perform the DFT and find the requested 
amplitudes Ej[0], and £y[l]. Of course, a si mple arithme t ic mea n 
can be used instead of the DFT (as we did in lBreiteretalJ ( l2009l) ). 
but the complete spectrum is required to compute the conductivity 
terms, as described in the next section. 



However prudent .seems the approach described by Scheeres (2007), de- 
tecting the longitudes of Sun path intersections with facets edges, it may 
involve another kind of subtleties responsible for the differences between 
Scheeres et al. (2008) and Breiter et al. (2009) concerning the influence of 
shadowing on the YORP eifect for 25 143 Itokawa. 



4 CONDUCTIVITY TERMS 
4.1 Plane-parallel model 

The surface temperature gradient, required by the conductivity 
complement, is obtained by solving the heat diffusion equation 
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dT 

V-KVT =pcp — 



dt 



(30) 



: F2T. 



(39) 



where p is the density and Cp is the specific heat capacity of the 
object. If we assume that conductivity K is independent on temper- 
ature and has the same value in the entire volume of the body, we 
reduce Eq. ([30} to the form 



(31) 



(32) 



dT 

kAT = — , 
dt 

where the thermal diffusivity k, defined as 
K 

will be assumed constant, leading to the homogeneous body ther- 
mal model. 

The plane-parallel model (PPM) results from two simplify- 
ing assumptions: i) the penetration depth of the heat wave is small 
compared with the radius of curvature for all fragments of the body 
surface, and ii) there is no heat exchange in the direction perpen- 
dicular to the surface normal. Both the assumptions are plausible 
for large objects with a low conductivity and a smooth, preferably 
convex surface. In PPM we introduce the depth variable f whose 
values increase from ^ = on the surface to higher positive values 
inside the body. The basic equation of the homogenous body PPM 
is a reduced form of ( 13 U 

d'^T _ dT 

with nonlinear Robin boundary conditions 
£t 7^(0) - r'(0) - ( 1 - A) £■ = 0, 



(33) 



(34) 



on the surface, and Neumann boundary condition in the limit of 
infinite depth 



lim T'iO = 0. 
In both cases we use 

and Eq. i34\ results from the energy balance ^ with 
hWT = -T'. 

Accordingly, Eq. il5i can be replaced by 



Qj 



KT'j. 



(35) 



(36) 



(37) 



(38) 



Instead of initial conditions at some specified epoch t, we 
impose the quasi-periodicity condition, requiring that all transient 
terms have been damped after sufficient relaxation time. This con- 
dition is most easily imposed by assuming from the beginning that 
T is replaced by its DPT with respect to rotation phase Q. and mean 
anomaly C. 

Since the assumptions of PPM exclude the heat transfer be- 
tween adjacent triangles (and their associated volumes), we may 
consider each body fragment separately, so the index referring to a 
particular face (like j in Eq. l l38t ) will be omitted in the following 
discussion. 



4.2 Helmholtz equation and its solution 

4.2.1 General case 

Let us consider the DFT of temperature 



Resorting to the associated trigonometric polynomial substi- 
tuted into Fourier equation ( I33t . we find that the DFT coefficients, 
as functions of depth obey a system of decoupled ID Helmholtz 
equations 



[f[p])"-ifipf[p]=0, 



p = 0,...,N^-l, 



(40) 



where A' is the angles sampling density, and parameters Pp depend 
on orbital mean motion v and rotation rate a>; if p = j + kN, then 



(41) 



where is defined in Appendix lAl 

Using a formal analogy with harmonic oscillator (with a com- 
plex frequency), and imposing the Neumann condition i35l we ob- 
tain a solution, depending on one arbitrary constant C^, in a form 



T[p] = Cp exp 



-(l + sgnC8,).) M( 



(42) 



In principle, Cp should now be determined from the second bound- 
ary condition, but for the further treatment we need only the loga- 
rithmic derivative 

^'' = 5t=-(1"^^S"M#' (43) 
T[p] ' / 

which occurs to be independent on <f and allows to express the 
derivative T' in terms of T. 



4.2.2 Null frequency 

The general solution J42t is not valid for pp = 0, when Eq. l l40t 
degenerates into 

[tip])" = 0. (44) 

It happens when p = 0, i.e. for the mean value of temperature T. 

The solution of l l44t is a linear function of f , but matching it 
with the boundary condition l |35l l we find that the null frequency 
solution is T[0] = const, hence 



70 = 0. 



(45) 



The fact that 70 = 0, has significant implications for the YORP in- 
fluence on (J). 



4.3 Boundary conditions 

4.3.1 Newton-Raphson setup 

Knowing the ratios jp, we can find the spectrum Q from the bound- 
ary conditions l l34t . Consider the vector of sampled temperature 
values T at the centroid of a given triangular face. We will desig- 
nate by T" the vector of the n-\h powers of T, i.e. 



T"[p] = (T[p])", 



p = 0,...,N^-l. 



(46) 



Then, using the DFT formalism from Appendix 1X1 the boundary 
conditions l |34| l can be written in the vector form 



etff-r4-e = (l-A)£. 

Using Eqs. l|A8ll, {MJ. CHl and lUll, we find that 

tUn-^f;bf2)t=—e, 



(47) 



(48) 
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where B is an x diagonal matrix witli 
£tcr ■ 



B 



(49) 



Tiie main difficulty in dealing with the energy balance equa- 
tion (|48} is its nonlinearity, requiring the use of some approximate 
methods. Resorting to the Newton-Raphson method, we can estab- 
lish an iterative scheme 

(D('"> + C)r*'"+"=G<'">, (50) 
where D is a diagonal matrix with 

<'=4(r3[p]f\ (51) 

C is a normal, block circulant matrix 

C = A'-2f*BF2, (52) 
with all eigenvalues Ap = Bpp having non-negative real parts, and 

(53) 



Q(m) ^ 1 ^ ^ ^ 2 □<'") 
6tCr 4 



In principle, starting from any reasonable approximation T^^\ we 
can solve the linear system J50b . approaching a sufficiently accurate 
T with a quadratic convergence. Then the spectrum Q, required in 
I l22l24t . follows from 



(54) 



efficiently executed by one call of the fast Fourier transform (FFT) 
routine. 

Regretfully, the left-hand side of Eq. ( 150b contains a dense, 
^2 ^ j^l jjj^fjjx carmot be directly inverted by low cost algo- 
rithms. This is quite frustrating, because the inversion of the diag- 
onal matrix D is trivial, while inverting FjB F2 alone is easily done 
by the FFT. But before we show the way to solve this problem, an 
important property of the nonlinear system | |48I > is worth stating. 



4.3.2 Conductivity has no influence on the rotation period in the 
PPM 



According to Sect. 14.2.21 and Eq. ( |49t , the element Sqo = 0. As a 
consequence, the first row of the matrix FjB contains only zeros, 
hence 



e[0] = 0. 



(55) 



Thus, returning to Eq. ( I22l ( we conclude that, as far as the plane- 
parallel model of a homogeneous body is concerned, the conduc- 
tivity complement has no effect on the mean value of the YORP 
torque component responsible for io, i.e. 



3 c 



(56) 



is determined by the Rubincam part alone. 

A similar observation was reported in previous works , 
alffiough each time with di fferent assumpt ions. iMvsenI j2008l). 
iNesvomv & Vokrouhhckvl ( l2008al) , and iBreiter & Michalskal 
( l2008h found it assuming the infinite radius of a homogeneous 
object, but they made additional assumptions of linearized temper- 
ature vari ations and pseudo-convex shad owing model. Numerical 
results of ICapek & Vokrouhlickvl ( l2004h . using the assumptions 
similar to these of the present paper, were nonconclusive: some 
objects seemed to have 6j independent on conductivity, but some 



(like 6489 Golevka) were exceptions from this rule0 The argu- 
ments base d on the spectrum of derivative support our earlier 
conjecture ( iBreiter et al]|2009h that the apparent dependence on K 
is definitely due to numerical errors - most probably a too short 
relaxation time and/or inac curate discretisation in the time stepping 
finite difference scheme o f lCapek & Vokrouhlickvl ( [2004h . 

The first significa nt dependence of u) on K was announced in 
the analytical model of IBreiter et al. 1(^010') which allowed for a fi- 
nite body radius and used a 3D heat diffusion equation in spherical 
coordinates, although - as usually in analytical models - with many 
additional simplifications. Thus, a central question is which of the 
two factors generates the dependence on conductivity. Appendix ICl 
presents the extension of the PPM to the ID model with finite ra- 
dius; even in this generalization 70 = 0, hence we can state, that 
the necessary condition for the dependence of the YORP effect in 
spin rate on conductivity is the heat exchange perpendicular to the 
surface normal, i.e. a 3D heat diffusion model. 



4.3.3 HN iterations 

In order to so l ve Eg . ( |50t . we took the approach based upon the idea 
of lHo&Nd ( 12003) who considered circulant-plus-diagonal sys- 
tems with imaginary diagonal part (iD-l- C). Unfortunately, major 
part of the proofs given by Ho and Ng relies on the skew-Hermitian 
property of ID, so we adopted their method to our (D + C) system 
faute de mieux, hoping that HN iteration j3 will work anyway. 

According to the HN algorithm, at each Newton step i50\ of 
the 'outer iterations', one should introduce 'inner iterations' 

(rl + C)y<** = (Tl-D("'>)r('«*+G<™', (57) 

(rl + D<"'>)r('"+l'*^ = (tI-C)F(*>+G<'"\ (58) 

where t > is some arbitrary real parameter, and Y e M.^'^ is an 
auxiliary vector. 

Concatenating Eqs. ( I57t and l l58t , one can see that the con- 
vergence of this process depends on the spectral radius p(M) of the 
matrix 



M = (tI + D)"' (tI-C) (rl + C)"' (rl-D), 
(superscripts (m) omitted) which is bounded by 

\T-Dpp\ \T-Bpp\ 

p(M) < max ■ — : max 



(59) 



(60) 



P \T + Dpp\ p \T + Bpp\ 

0, or %(Bpp) = \5(Bpp)\ > 0, we con- 



Knowing that either Bpp 
elude 



\r-Bpp\ 
max 

p \T+Bpp\ 

hence 



1, 



p(M) < max 



\t-D„ 



P \T + Dpp\ 



(61) 



(62) 



Thus the spectral radius is less than I, provided the diagonal of 
D does not contain zero or negative values of 4T^. It means ffiat 
HN iterations will converge faster at higher conductivity values, 
when the minimum temperature does not drop significantly during 
the night. On the other hand, Eq. J62t suggests a safe and nearly 



2 ICapek & Vokrouhlickvl f2004 write about 'a near independence' on K. 
^ An acronym equally matching the authors' names and the Hermitian- 



plus-normal nature of the system. 
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optimal choice of r as the geometric mean of the maximum and 
minimum diagonal entries of D 



1"= VAnaxAi 



(63) 



This rule obviously fails for D^[„ = 0. But, what is worse, the shad- 
owing effects lead to discontinuities in the insolation, causing the 
so called ringing artifacts - often with negative values of temper- 
ature. From the point of view of upper bounds l l62t . the iterations 
should diverge in such cases, but the algorithm occurs to be un- 
expectedly robust and often converges in spite of T < 0, although 
once the temperature drops below 0, a number of wild and chaotic 
jumps can be observed before the residues resume their decreasing 
path. After a number of trials we have finally adopted a practical 
rule of thumb 



■ = max(VD max 



(64) 



handling the negative Z?,nin case, and protecting t from taking ex- 
cessively small values. 



t('">-S„„ , 



(72) 



pp 



and then we obtain the m-th approximation of surface temperature 



T-(m) + 

' ^ PP 



(73) 



Each step of this process requires two calls of the FFT pro- 
cedures, one direct (X — » X) and one inverse (W — > W), as well 
as few loops with n'^ complex products, which is probably not far 
from the optimum computational cost. 



4.4 First guess and accuracy 

The fundamental question accompanying each iteration process is 
how to start and when to stop. It looks reasonable to assume the 
starting value T^^' by setting K = Q 'in the original, nonlinear bound- 
ary conditions ( |48l >, which leads the choice between 



4.3.4 Quasi-Newton method 

When using combined inner-outer iterations schemes, one always 
faces a problem when to terminate the inner loop before the im- 
provements become nonsignificant from the point of view of the 
current outer iteration. At this point we trade efficiency for sim- 
plicity and retain only one inner HN step, obtaining the final quasi- 
Newton scheme with two substeps 



(65) 
(66) 



(r^™^ I + C) F = (r^'"^ I - D<'"') T*'"' + G^"'\ 

(j(m) I ^ Q(m) ^ 1 ) ^ (^^(m) | _ c) F -I- G*'"' , 

where t('«', D('"\ and G<"'' are recomputed at each step m (but not 
between the substeps J65b and J66b). 

In practical terms, solving the equations of the quasi-Newton 
method is quite simple and requires the storage of only few ID 
arrays with elements. The matrix- vector product in the right- 
hand side of Eq. ( 165 1 is obviously executed in a single loop, 
generating some vector X 6 R'^ according to 

X[p] = (t("') - d(,7,') r('«) [p] + GC") [p], (67) 

where p = 0,...,N^ - I . We compute this vector and find its DFT 
X according to the definition iA6l . In order to solve the system 
(t("') I -I- C) F = we note that 

(t*'"' I + C) F = (F* [n-^ l) F2 + C) F, (68) 

so, substituting Eqs. l |52||67|[68l l into ([65} we obtain 

F;(t<'"'| + B)f = Z, (69) 

where F = A'"^F2F is the DFT of F. Thus, the first substep is com- 
pleted by defining, but not yet evaluating, the transform F 



Y[p] ■■ 



X[p] 



r('") + Br, 



(70) 



Solving Eq. \66\ we use a similar approach: first, the product in the 
right-hand side is expressed as 



(t(™'|-C)f=F5(t('"'|-B) F. 



(71) 



It means, that we have to compute the inverse DFT W = FjW, 
where 



y(0) ^ 



£t(T 

or, apparently simplistic. 



rW[p]=(i-A) 



6tcr / 



-To, 



P = 0, 



(74) 



(75) 



Choosing a constant T*^"' according to ( I75t may seem too crude, 
since it means that iterations will have to reconstruct all periodic 
terms with leading amplitudes - in the worst case of low conductiv- 
ity - comparable in magnitude to the mean value. But the practice 
shows a superiority of ( I75t over ( I74t . Building the amplitudes up 
from zero is numerically more stable than decreasing their values 
from the state, when the temperature determined by ( 174b takes zero 
values. This fact can be explained from a number of points of view, 
using both physical and mathematical arguments. Focusing on the 
latter, note that according to the estimates given in Sect. 14.3.31 the 
spectral radius p(M) equals 1 when any of T^^^[p] - 0. Moreover, 
the approximation ( 174 1 is a continuous, but not smooth function of 
t and fi, which significantly degrades the numerical quality of the 
DFT of its derivative with respect to these angles. 

Let us write explicitly the values of T resulting from the 
quasi-Newton iterations when T*''' is given by Eq. ( I75t . In this case, 
/"'"^[p] = for all p 7^ 0, the mean value is 7'**'^[0] = Tq, and diag- 
onal matrix D(°> = 4r^ I, hence t*'^' = 4r^. Using 



GW = iZii£ + 37-3 7(0), 
etO" 

we obtain from J65b and J66t . left multiplied by F2, 



r = \ — 



p=l,...,N^-l. 



(76) 

(77) 
(78) 



'^"PP 



Remarkably, the same result can be obtained even easier from the 
original Newton-Raphson system ( I50t . Equations ( I77t and ( I78t de- 
fine the linear thermal model - a standard tool in analytical YORP 
theories. Of course, the direct application of Eqs. ( |77t and ( |78l l, 
followed by one FFT call to obtain T^'' is much cheaper than per- 
forming the first iteration of ( 165 1 and ( |66l > in its regular form. For 
this reason we actually start iterations from m = I, obtaining the 
means to simulate the results of linear model as an extra profit. 
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Table 1. Parameters of test objects 



1998 KY26 









1998 KY26 


6489 Golevka 


albedo 


A 




0.1 


0.1 


emissivity 


ft 




0.9 


0.9 


density 


P 


kgm ' 


2800 


2700 


specific heat 


Cp 




680 


680 


mom. inertia 


C 


kgm^ 


1.9346x10' 


7.3420 xlO'5 


rotation period 


P 


h 


0.1748 


6.0264 


orbit semi-axis 


a 


au 


1.232 


2.5 


eccentricity 


e 




0.2 


0.6 


approx. diameter 




m 


~ 30 


-530 



The iterations cycle has to be stopped when a sufficient accu- 
racy is attained. The stopping criterion should be chosen carefully. 
The simplest one is to observe the differences between subsequent 



values of 7'*'"'[1] and exit when |7'''"'[1] - 



f^"' '^[1]|<(5. Butone 



has to be careful, because when the convergence is slow (a typical 
situation at low conductivities), such a difference carries no infor- 
mation about the accuracy of 1]. Fortunately, we have also an 
objective criterion eicrlr^'^ = (1 -A)£'[0], independent on the con- 
vergence rate. And since the mean value of T"* accumulates also 
the errors of all periodic terms of T, we exit the iterations when, for 
a specified temperature tolerance 6, two conditions are simultane- 
ously satisfied; 

-3 llr^(m)\^\ (1- 



-A) - 

■Em 



<6. 



(79) 



5 TEST RUNS 

5.1 Test bodies and accuracy requirements 

Two asteroids were chosen as test bodies for our numerical sim- 
ulations: 1998 KY26 with a relatively regular shape, and 6489 
Golevka, whose large scale concavities and sharp corners make 
it a good benchmark for the YORP models. Physical parameters 
adopted for the simulation are given in Tab. [T] Radar shape models 
of both object!0 (4092 triangular faces) were reduced to the center 
of mass and principal axes system assuming a constant density. 

In most of previous YORP models either the orbits were sim- 
ply assumed circular, or the YORP effect computed on a circular 
orbit was multiplied by a conversion factor 



(80) 



From theoretical standpoint, the latter procedure can be justified ex- 
act in the Rubincam's approximation or in linear thermal models, 
where rotation and orbital motion effects are separable, but there 
are no reasons to assert it in general. The factor concerns all 
terms proportional to the average power flux (i.e. the ones with 
(y*)) but not tho se depending on the first power of temperature 
( lRubincam|[200 4'). 

In all computations we have adopted a rule that no error 
bars should be required in YORP plots. Various levels of sam- 
pling A' and tolerance 6 had been tried until a difference from 
the results with sampling 2N and tolerance 6/10 became compa- 
rable to the plot line thickness. Finally, for 1998 KY26 we used 



* Downloaded from the site http : //www . psi . edu/pds/asteroid/ file 
version EAR^^.DDR.RADARSHAPEJVlODELS.V2.0.zip 
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Figure 2. YORP effect in rotation rate for 1998 KY26 (top) and 6489 
Golevka (bottom). Joined dots mark the values computed with 2° spacing 
in obliquity. Dashed line represents the pseudo-convex approximation. Note 
the difi'erence in units between the top and bottom panels. 



N = 256, (5 = 10 , whereas Golevka, as expected, was more chal- 
lenging and required A' = 512, 6 = 10 , or even lO"*", depend- 
ing on conductivity. For the pseudo-convex shadowing model the 
requirements were less severe and A' could be two times smaller. 
Formal accuracy of the results presented in next sections is the fol- 
lowing: 1998 KY26 ~ 3 X 10"^ rads^' My"' for the rotation rate, 
and 3x lO"'* rads^'My-i for the attitude YORP; 6489 Golevka 
- 0.03 radd^'My"' for w, and 0.5 radd^'My"' for the attitude 
effect. 



5.2 YORP effect in rotation rate 



As we demonstrated in Sect. 14.3.21 and Appendices |B] and [C] all 
kinds of ID thermal models lead to the same YORP effect in ro- 
tation rate, equivalent with the Rubincam approximation A" = 0. 
Figure |2] shows the values of doubly averaged (M3)C"', where 
M3 = M-e^. According to Eq. l ll6t . these values are equal to angu- 
lar acceleration di. The dots in Fig. |2] are placed for actually com- 
puted values of co, and they form curves tha t fairly well agree with 
the results o f IVokrouhfickv & Capekl ( |2003) . provided the factor 
is used and the differences in C and are accounted for. 

The pseudo-convex approximation looks decent for a regu- 
lar object like 1998 KY26, but it fails completely for irregularly 
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shaped Golevka. It is worth noting, that the influence of shadow- 
ing for 1998 KY26 amounts in principle to a vertical translation of 
the curve; a similar (although more prominent) phenomenon was 
observed for 25143 Itokawa (.Breiter et al.,,2009.) . 



5.3 YORP in attitude: seasonal effect revealed 

The YORP effect in attitude is usually described in terms of 
(M-ei) = (Ml), and (M-e2> = (Af2>- According to Sect.|2j2l the 
mean value of the drift in obliquity is given by oje = (Mi)C"', 
whereas (M2)C"' = tane^fl-o)) is responsible for the mean pre- 
cession component of the effect. Figures|3]and|4]demonstrate the in- 
fluence of conductivity on these two attitude components within the 
nonlinear ID thermal model. Ex cept for the Rubincam case, wh ere 
a good agreement of (Mi) with IVokrouhlickv & CapeS ( l2002h is 
observed, the results may look surprising, not to say ridiculous, at 
the first glance. Isn't it absurd to have e > for e = 180° ? Why are 
the present curves so different from all previously reported plots ? 

The first question is relatively easy to resolve, since it is re- 
lated to the classical problem of polar coordinates singularity close 
to the origin, where a wrong parametrization may contradict physi- 
cal facts. Nonzero mean values of (Mi) and (M2) at sine = merely 
indicate that the orientation of the spin axis normal to the orbital 
plane is not an equilibrium. A proper treatment of passing through 
this state requires a formulation in terms of the spi n vector and 
torque Cartesian coordinates (e.g. lSreiter et al.l j2005l) ). 

As for the second question, we have to note that the previ- 
ous theories of YORP with nonzero conductivity were mostly lin- 
ear, approximating T"* by + 4T^ Ti , with a constant Tq and a 
purely period ic Ti . The on ly exception from this rule is the numer- 
ical model of ICapek & V okrouhlickv (2004), but there the authors 
present only the results for circular orbits. They did compute the 
values for e as well, but only for single, specific e and oj^ pairs 
of actual objects and no plots covering the whole range of obliqui- 
ties have been published as yet. Figure [5] shows that linear approx- 
imation generated by our model (top) and nonlinear results with 
e = (bottom) behave exactly like in previous publications (except 
for a more complicated shap e resulting from a better sampli ng than 
the 9 points interpolation of ICapek & Vokrouhlickv! ( |2004|) ). Even 
a weak asymmetry o f the obliquity YORP curve with respect to 
e = 90° agrees with dCapek & Vokrouhlickv! I2OO4I) . These results 
imply that the shape of curves in Figs.|3]and!4]is due to nonlinear 
coupling between daily and seasonal waves, and the effect must be 
due to the variation of heliocentric distance, because the tempera- 
ture variations due to change of seasons are present also in circular 
motion where nothing unusual happens. 

More light can be shed on the problem when changing the ar- 
gument of perihelion 0)^, which was set to in all previous plots. 
Figure |6| presents the attitude YORP effect for Golevka (e = 0.6, 
K = 10"^ Wm"' K"') with four different arguments of perihelion 
(1)0 values (0, 90, 180, and 270 degrees). Figure |7| compares the 
arithmetic mean of the four values with the results obtained for the 
circular orbit and re-scaled to e = 0.6. On the other hand. Fig. |8| 
shows the dependence of the attitude YORP effect on the argu- 
ment of perihelion (sampled by 4°) when we fix the obliquity of 
Golevka at £ = 30° . The dependence is almost (but not exactly) si- 
nusoidal and the amplitude depends on eccentricity, although the 
dependence does not seem to obey a simple power law. 



5.4 The two layers model 

Introducing a more advanced model with a monolithic core and a 
regolith layer, described in Appendix !B] has no influence on the 
YORP effect in rotation rate. Thus, the results shown in Figs. |9] 
and 1 101 concern only the effect in attitude. In the test runs we have 
considered the values of density for both objects given in Tab. [T] 
as the bulk densities serving to compute the moments of iner- 
tia, but they are no longer used to compute the thermal diffusiv- 
ity. Instead, we have adopted the following regolith parameters: 
K = O.OI Wm''K' ', cp = 760 Jkg^' K^', and p = 1660 kgm^^ 
( iRumpf et alj !2008h . Accordingly, the thermal diffusivity of re- 
golith layer is /c 7.93 x 10"^ rc? s"' and the physical properties of 
the core are specified by the ratio w defined by Eq. l lB3t . We have 
assumed the value of w = 0. 1 as a presumably realistic estimate. 

Fixing a randomly chosen obliquity e = 45°, we used our 
model to compute the YORP effect in attitude for various regolith 
depths h. As it might be expected, there is a gradual transition be- 
tween the thin and thick regolith cover results, and the curves in 
Figs. |9] and [Toj are practically flat when prolonged towards higher 
or lower h values. However, the transition is not monotonic, resem- 
bling a superposition of a logistic curve with damped oscillations. 
This effect is understandable, observing that h factors both the real 
and the imaginary parts of the exponential in Eqs. l |B8|B9t . A sim- 
ilar pa ttern was present in the Yarkovsky force model of !CapelJ 
( !2008i Fig. A.8) - the only analogue that we can refer to. 

The characteristic order of magnitude for the depth h deter- 
mining the transition from the thick to thin regolith case is the skin 
depth: a function of thermal diffusivity and insolation frequency 
[Lagerros!! 19961) . However, there are two different principal skin 
depths in our model: rotational involving ai, and orbital Zq, in- 
volving the mean motion v: 

Vertical lines in Figs.!9land|10|mark these two parameters (dashed 
for /r and dotted for /o), indicating that rotational skin depth (much 
smaller than /q) is the only important quantity. However, a signifi- 
cant deviation from the thin regolith mode occurs already at the val- 
ues of h below 0. 1 k or even O.OI k. With ~ 1 mm for 1998 KY26, 
and Ij X 5 mm for Golevka, we can observe that the strongest de- 
pendence of the attitude YORP effect on h is observed when the 
regolith thickness is in the range of O.I -r 10 mm, wh ich is quite 
similar to the results of ! Vokrouhlickv & Broz ! ( !l999!) concerning 
the Yarkovsky effect. 



6 CONCLUSIONS 

Thanks to the application of Fourier transform, the algo- 
rithm presented in this paper is more efficient and more ac- 
curate (although less genera l) th an its equival ent described by 
!Capek & Vokrouhlickv! ( !2004!) and !Capek! ( !2008!) . Abandoning the 
finite difference approach in favor of using exact solutions of the 
Helmholtz equation helped us to demonstrate, that the YORP ef- 
fect in rotation period is the same in the Rubincam's approximation 
(K = 0) and in various ID thermal models. As long as one neglects 
the heat transfer between adjacent surface elements, the values of 
(1) do not depend on conductivity, regardless of the body size and 
radial homogeneity assumptions. From the point of view of obser- 
vational detection of YORP, always based upon di, this is a nice 
conclusion; not only because the Rubincam's model is easier to 
compute, but also it requires less physical parameters to be known. 
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Figure 3. YORP effect in obliquity (left) and precession (right) for 1998 KY26 with e = 0.2 and oj„ = 0. Joined dots mark the actually computed values in 
Rubincam's approximation (K = 0). Dashed, dash-dotted and dotted curves refer to the values from nonUnear ID model at K = 0.001, 0.1, and 10 Wm'K"', 
respectively. 
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Figure 4. Same as Fig.[3]for 6489 Golevka with e = 0.6 and a)a = 0. 



since then - at least for the Lambertian scattering and emission - 
emissivity and albedo values do not matter. As a matter of fact, 
the conclusion can be also given a straightforward physical expla- 
nation. If the heat conduction is restricted to the direction normal 
to the surface, a nonzero mean value of the the temperature normal 
gradient Q should imply systematic heating or cooling of asteroid's 
interior. Hence, the property that (Q) = follows directly from the 
request of the energetic equilibrium state with transient te rms re- 
laxed . However, according to the analytical model of iBreiter et al.l 
( I2OIQ) . the YORP torque in rotation period for smaller bodies with 
a 3D thermal model may differ from the Rubincam's approxima- 
tion because of the heat flow between adjacent surface pathes that 
may receive a different mean power flux. 

Even if the YORP effect occurred to be insensitive to the 
transverse heat conduction, the ID models considered in this pa- 
per should not be applied to objects whose diameter is small when 
compared with a skin depth. This restriction, explicitly stipulated 
in Sect. |4.1| can be physically explained as follows. Consider a bar 
passing through the centre of a body O and intersecting the surface 
in two antipodal areas S i and 52- The ID models consider it as two 
disjoint slabs with the absence of heat conduction at O imposed as 
a boundary condition. In these circumstances, even if the conduc- 
tivity is very high, there is no possibility to transfer the heat from 



the sunlit S 1 to the dark 5 2 in order to reach a smoother surface 
temperature distribution and reduce the YORP strength. A possible 
improvement of ID models might be based upon considering the 
set of antipodal bars without the central cut; yet, in our opinion, a 
future investment in a complete 3D model is more needed. 

The YORP effect in attitude is not directly observable, but 
still important for the simulations of long-term spin axis dynam- 
ics. The most prominent example is the analysis of the Slivan states 
in Koronis family ( Vokrouhlickv et al. 2003), considered the first, 
indirect proof of the YORP effect existence and significance. Our 
results indicate that for elliptic orbits there exists a phenomenon 
that may be called a seasonal YORP effect in attitude by analogy 
with the seasonal Yarkovsky effect in orbital motion dRubincaml 
1995; Vokrouhlickv & Farinella 1999). The seasonal effect did not 
appear in earlier works based upon linearized thermal models, 
which led to a hasty rule that the influence of orbital eccentric- 
ity amounts merely to a multiplicative factor from Eq. l lSOt . We 
confirm the validity of this rule for the rotation period YORP, but 
not for the attitude. The effec t passed unnoticed in the model of 
ICapek & Vokrouhlickv! ilOOA which can probably be explained 
by its high computational time demands that discouraged exper- 
iments with various argument of perihelion values. The seasonal 
YORP in attitude deserves a closer inspection within a nonlinear 
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Figure 6. YORP effect in obliquity (left) and precession (right) for 6489 Golevka on an eccentric orbit (e = 0.6) with = (solid), 90° (dotted), 180° 
(dot-dashed), and 270° (dashed). Conductivity K = 10"^ Wm"' K"' . 



analytical model (even with a crude insolation model) that might 
help to explain its physical meaning. It is quite possible that there 
exists some relation between the seasonal Y ORP and the Sev er- 
smith psychroterms mechanism discovered bv lRubincanil ( l2004h 

As a final remark, let us observe that the presented model can 
be easily adapted to compu t e the Yarkovsky effect in orbital motion, 
like in the paper of lMvseiJ JioOSf) . 
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Figure 7. Arithmetic mean of tlie four curves from Fig.|6](solid line) and the results for a circular orbit multiplied by 1.25 (dots). 
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Figure 8. YORP effect in obliquity (left) and precession (right) for 6489 Golevka on an eccentric orbit (e : 
10"^^ Wm"' K"' as a function of argument of perihelion ai^. 



0.6) with £ = 30° and conductivity K - 
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Figure 9. YORP effect in obliquity (left) and precession (right) for 1998 KY26 at £ = 45° as a function of regolith depth h. Rotational and orbital thermal 
penetration depths are indicated by vertical lines (dashed and dotted, respectively). 
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Figure 10. Same as Fig.|9]for 6489 Golevka. 



APPENDIX A: DISCRETE FOURIER TRANSFORM 

Our implementation of the algorithms presented in this paper relies 
on the extensive use of the FFTW library (v. 3.2.2) developed by 
lFrigo& Johnsoi] ( l2005l) . The following formulae will use sign con- 
ventions, normalization factors and ID storage of 2D matrices (in- 
cluding the numbering of elements from 0) adhering to the FFTW. 
The size N Fourier matrix F is defined in terms of powers of 

wy = e-'>2^/^, i = Q,...,N-\, (Al) 



,N-\ 



(A2) 



For the 2D discrete Fourier transform (DFT) of the size NxN,we 
define the matrix F2 using the Kronecker tensor product 

F2 = F(8iF, (A3) 

with a resulting block structure of the x A'^ matrix 

^0 ^ 

(A4) 



, ,0 P 



,N-\ 



F j 



Let us consider a function of two angles j((0, 1/'). Sampling 
u on a square grid of 0j = jlnjN, and tp^. = kin IN, where j,k = 
0, . . . , W - 1 , we create a vector u e WJ^'^ , whose p-th element is 



u[p] = u{<pj,\l/k), p = jN + k. 



(A5) 



The direct discrete Fourier transform (DFT) of u is the vector 
u 6 M.^ resulting from the matrix-vector product 



u = — ^ Fou. 



(A6) 



The inverse DFT is provided by the complex conjugate Fj with the 
property 



F*F2 = F2F; = A'2|, 
so that 

Ft « = — t-F*Fm = M, 

2 N2 



(A7) 



(A8) 



explaining the necessity of the N^- factor in Eq. l |A6l >. 

We can consider DFT as the coelficients of a trigonometric 
polynomial 



A'' N' 

Z Z N' = lN/2i, 

j=-N' k=-N' 



(A9) 
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(AlO) 



where L J is the "floor" rounding down operator. Introducing 

7 , , / 9 for q^l^Nl 
^^^^^ = 1 for q>llNi, 

we can identify 

u[qiN + q2] = Ujk, (All) 

j = ZNiqi), (A12) 

k = ZNiqi), (A13) 

with the indices q\,q2 = 0,...,N- 1. Strictly speaking, for even 
the Nyquist terms with |j| = N/2 or l^-] = N/2 require a special 
treatment and an extra factor 1/2 or 1/4, but their influence on the 
final solution is practically so marginal, that we do not pay attention 
to this problem. 



APPENDIX B: PPM WITH REGOLITH LAYER 

The thermal model presented in Sect. |4] can be easily extended to 
cover a case when an asteroid is treated as monolithic core covered 
by a regolith layer of thickness h. Let us assume conductivity K and 
thermal diffusivity k for the regolith layer < ^ < A, and Kc, for 
the core f h. Then, Eq. (|40} is replaced by two sets 



[hp})' -i/ipTip} 

[f,[p])' -iw^/3pf,[p} 



(Bl) 
(B2) 



where p = Q,...,N -I, and w is the ratio of thermal diffusivities 



(B3) 



introduced to use a single parameter /3p defined for the surface layer 
according to Eq. (l41tPl 

Solutions Ti^[p\ are subject to the Neumann condition at infin- 
ity, hence, similarly to ( I42t they read 



TAP^ = Cp exp 



-(l+sgnCe;,)i) 



w\/3p\ 



(B4) 



except for the special case 7'c[0] = Cq. For the regolith layer, how- 
ever, the asymptotic condition does not apply, so it takes a more 
general form T[p] 



Tip] 



Ap exp 



+Bp exp 




(B5) 



with the boundary condition l |48l ). As usually, the special case 
7'[0] = Aq applies. 

Both solutions should satisfy continuity requirements at( = h, 
i.e. additional Dirichlet conditions 



{fc[p]-f[p])^^^^=0, 
and Neumann conditions 

{K[p]-f'[p])^^^ = o. 

Our w is the same as of lVokrouhlickv & BrozNl999l) . 



(B6) 



(B7) 



Thanks to them, we can express Bp in terms of Ap alone, and so 
we obtain the logarithmic derivative that generalizes ( 143 1 for the 
surface temperature for p ^ 



yp 



f'lp] 

hp] 



-(l+sgnOSp)i)y 



w-l + (w+ l)Hp 



where 



Hp=e 



M (B8) 
■w + (w+l)Hp V 2 



(B9) 



and the special case with p = is 70 = 0. 

The formula JB8t uses no assumptions about the depth h, but 
if we postulate a thin regolith layer, the terms linear in h are simply 



- (1 + sgnC8p)i) wJ^+ i/3p (w2 - 1) h. 



(BIO) 



Recalling that w pp is the ratio of insolation frequency to the ther- 
mal diffusivity of the core, we see that jp is dominated by the prop- 
erties of the inner part of asteroid, although entering the matrix B it 
will be multiplied by the surface conductivity of the regolith layer 

a:. 

Going even further and neglecting h in the approximation 
l IBlOt . we obtain the simplest recipe for the regolith covered ob- 
jects: use the physical parameters of the core to compute yp, and 
the regolith parameters in the rest of the algorithm. Of course, it 
would be wise to verify the quality of this approximation by com- 
paring at least the values of 71 computed according to l IBlQI l with 
h = 0, with the ones obtained using a more exact formula. 

Concluding this section we emphasize that even when the 
physical properties of an object vary with the depth i^, our con- 
clusion about the independence of (w) on conductivity remain true. 



APPENDIX C: SPHERICAL SEGMENTS AS A ID MODEL 
WITH FINITE DEPTH 

Suppose that an object is starlike, i.e. each surface element can be 
connected with the centre of mass using a straight segment that 
does not intersect other surface elements. In these circumstances, 
we can consider the Laplacian operator in spherical coordinates, 
and the reduction to ID amounts to neglecting non-radial part of A, 
so that Fourier equation Olt becomes 



K_d_ 

dr 



. dT 

a7 



dT 

It' 



(CD 



where < r < An elementary substitution 

11= -T, (C2) 
R 

reduces Eq. iCll to the form similar to l l33t 
d^ii du 

so the associated Helmholz equation is the same as ( I40l l 

^^-i/3 u[p] = 0, p = 0,...,A?2_i (C4) 
dr^ 

Note, that at the surface, where r = R,we have ;/ = T, allowing us 
to use the algorithm of Sect.|4]directly, without even changing the 
symbols. But first we have to redefine the coefficients yp to account 
for new boundary conditions. 

Instead of asymptotic Neumann conditions at <^ — ^ 00, typical 
for the plane-parallel case, we now have the Dirichlet condition 
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Up(r = 0) = 0, satisfied automatically when T at the origin is finite. 
On the surface, where u = T, the energy balance l|4j holds true, but 
now the gradient of T, reduced to the radial derivative, is not the 
same as the derivative of u, because 



dT 
d7 



r=R 



du 
d7' 



u 



Thus, instead of l l38t we use 

(;( d;( \ 
K 
R dr 



(C5) 



(C6) 



where 6 is the cosine of the angle between the radius vector and 
the outward normal vector of the current surface element. As we 
see, there are two major differences in the conduction treatment 
between the PPM and the spherical segment approach: the depen- 
dence on the local radius R, and the deviation of the gradient from 
the normal direction. 

Similarly to the PPM model, we solve Eq. iC4\ as a harmonic 
oscillator with imaginary frequency, but this time a difi^erent bound- 
ary condition gives us at the surface r = R 

1 du[p] 



, coth 



u[p] dr 

and the ratios yp in matrix B should be replaced by 



J I 1 d«[p] \ 
\R fi[p] dr j' 



(C7) 



(C8) 



when p + 0. However, when p = 0, we still have 

70 = 0, (C9) 

so the YORP effect in rotation rate w remains insensitive to con- 
ductivity, similarly to the PPM case. 

Taking the outward normal parallel to the radius, i.e. fixing 
(5=1, and taking the limit at infinite R, when /jp — > 1, tp — > 0, we re- 
cover the plane-parallel model with 70 = 0, and remaining jp given 
by Eq. 



